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ABSTRACT 

We present Keck/NIRC2 and OSIRIS near-infrared imaging and spectroscopy of 2M044H-2301 
AabBab, a young (1-3 Myr) hierarchical quadruple system comprising a low-mass star, two brown 
dwarfs, and a planetary-mass companion in Taurus. All four components show spectroscopic signs 
of low surface gravity, and both 2M0441+2301 Aa and Ab possess Pa/3 emission indicating they 
each harbor accretion subdisks. Astrometry spanning 2008-2014 reveals orbital motion in both the 
Aab (0'.'23 separation) and Bab (0'.'095 separation) pairs, although the implied orbital periods of 
>300 years means dynamical masses will not be possible in the near future. The faintest component 
(2M0441+2301 Bb) has an angular iif-band shape, strong molecular absorption (VO, CO, H 2 O, and 
FeH), and shallow alkali lines, confirming its young age, late spectral type (LI ± 1), and low temper¬ 
ature (^1800 K). With individual masses of 200^50^ Vfjup, 35 ± 5 Mj^p, 19 ± 3 Mj^p, and 9.8 ± 
1.8 Mjup, 2M0441+2301 AabBab is the lowest-mass quadruple system known. Its hierarchical orbital 
architecture and mass ratios imply that it formed from the collapse and fragmentation of a molecular 
cloud core, demonstrating that planetary-mass companions can originate from a stellar-like pathway 
analogous to higher-mass quadruple star systems as first speculated by Todorov et al. More generally, 
cloud fragmentation may be an important formation pathway for the massive exoplanets that are now 
regularly being imaged on wide orbits. 

Subject headings: stars: low-mass, brown dwarfs — planetary systems — stars: individual (2MASS 
J04414565+2301580, 2MASS J04414489+2301513) 


1. INTRODUCTION 


Over the past two decades a variety of planet-finding 
techniques have uncovered an unexpectedly diverse array 
of orbital architectures among extrasolar planetary sys¬ 
tems. Giant planets with masses between 1 to 10 times 
that of Jupiter have been found spanning over five orders 
of magnitude in separation from their host stars— from 
hot Jupiters at 0.01 AU to planetary-mass companions 
with projected separations exceeding 1000 AU. Uncov¬ 
ering the origin and dynamical histories of these popu¬ 
lations has emerged as a fundamental goal of planetary 


Several formation routes have been proposed to ex¬ 
plain the disparate orbital properties of giant planets. In 
circumstellar disks, the growth of rocky cores and subse¬ 
quent accretion of gas is generally accepted as the domi¬ 
nant mode of giant p lanet formation at s mall separations 
within about 10 AU (jPollack et al.|1996|). On wide orbits 
of tens of AU to a few hundre d AU,"^s giants may as¬ 


semb le from pebble accretion (Lambrechts & Johansen 


2012) or fragment directly from gravitational instabili- 


ties m massive protoplanetar y disks dur i ng the earliest 


phases of disk ev olution (e.g., Boss||19^ 


Whitwor'th||2009 ). Collapsing molecular c 
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urally give rise to binary systems with companion masses 
that span low-mass stars, brown dwarfs (<75 Mj^p), 
and possibly opacity-limited fragments with masse s be¬ 
low the denterinm-bnrning limit of about 13 y.Tun (|Low| 


& Lynden-Bell||1976| |Bate|2009t |Brandt et al.||20i4|). All 
of these mechanisms theoretically can produce gas giants 
between 1 to 10 Jupiter masses, making it difficult to 
disentangle the origins of individual planetary systems. 
Observational tests of these models are further compli¬ 
cated by high extinction at young ages as well as orbital 
evolution through disk migration, dynamical scattering 
from other planets or stars, and secular effects like Kozai- 
Lidov oscillations, all of which gradually er ode the frag¬ 
ile fos silized signatures of planet formation (jDavies et al. 

2M0441+2301AabBab is a young (1-3 Myr) quadru¬ 
ple system in the Taurus star-forming region that of¬ 
fers valuable clues about the ability of cloud fragmen¬ 
tation to_^orm_objects__^low the deuterium-burning 
limit (jTodorov et al.|2010). Indeed, following its dis¬ 
covery Daied"'drrimagmg data from th e Hubble Space 
Telesc ope /WFFC2 and Gemini-N/NIRI, |Todorov et al. 
(|2010|) suggest this system is most consistent with for- 
mation from a fragmented cloud core because of its 
extremely young age and low mass ratios. This sys¬ 
tem consists of two close binaries separated by 1273, or 
about 1800 AU in projection at the characteristic dis¬ 
tance of Taurus (145 ± 15 parsecs). The more mas¬ 
sive pair, 2MASS J04414565+2301580 Aab (hereinafter 
2M0441+2301 Aab), is composed of a low-mass star 
with an optical spectral type of M4.5 (the “Aa” compo¬ 
nent) orbited by a brown dwarf (“Ab”) at 0723 (33 AU). 
2MASS J04414489+2301513 Bab is composed of an M8.5 
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Fig. 1.— Point spread function (PSF) decomposition for extracting spectra from our Keck/OSIRIS observations. Left: Examples of 
median-combined cubes for 2M0441-h2301 Aab in Jbb (top row, blue), Hbb (middle row, green), and Kbb (bottom row, red) filters. The 
left panels show the data, the middle panels show our model fits composed of three bivariate Gaussians for each component, and the 
right panels show the residuals of the fits. The images are displayed with a logarithmic stretch. The median rms of the residuals for all 
observations of the Aab pair are 0.13, 0.23, and 0.31 Data Numbers (DN) for Jbb, Hbb, and Kbb, respectively. For comparison, the average 
peak values of the median-collapsed cubes are 10.1, 23.8, and 24.8 DN. Right: Same as the left panels but for 2M0441-h2301 Bab. The rms 
values for the Bab pair are 0.0056, 0.0097, and 0.012 DN for Jbb, Hbb, and Kbb. The average peak values are 0.14, 0.27, and 0.51 DN. 


brown dwarf (“Ba”) with a fainter companion (“Bb”) at 
O'/l (15 AU) whose apparent magnitude implies its mass 
is only about 10 times that of Jupiter. The hierarchi¬ 
cal arrangement of orbits between the compact binaries 
and much larger binary-binary pair suggests the system 
is dynamically stabl e. Both pairs are w ell-established 
members of Taurus dLnhman e t al.||2QQ9), have lit tle or 
no reddening (Ay ^ 0.2 mag; |Bulger et al. 2014), and 
are confirmed to be physical binaries rather tha n chance 
alignments based on multi-epoch astrometry (Todorov 
et al.|2014 ). An excess of mid-infrared emission was pre- 
viously identified from both unresolved pairs, indicating 
that at least one compon ent of each subsystem harbors a 
circu m-(sub)stellar disk (Bulger et al.|2014 Adame et al 
^Mq| ). Spectroscopic validation is essential to confirm 


the low temperatures, luminosities, and masses of can¬ 
didate planets at these very young ages when edge-on 
disks around brown dwarfs or low-mass stars can cause 
extinction an d mimic the photometric properties of faint 
protoplanet s ( Kraus et al.||2Q14 Bowler et al.||2Q14 Wu 


et al. 


2015). 


Here^ we present 1.2-2.4 /rm spatially resolved 
imaging and spectroscopy of the four components of 
2M0441+2301 AabBab. Our observations are described 
in Section spectral classification, analysis of emission 
lines, and atmospheric model fits are detailed in Section 
1^ and we discuss the implications of our spectra in the 
context of planet formation models in Section]^ 

2. OBSERVATIONS 


We obtained moderate-resolution (A/A A ^ 3800) in¬ 
tegral field spectroscopy of 2M0441+2301 Aab and Bab 
with the OH-Suppressing In frared Imaging Spectrograph 
(OSIRIS; [Larkin et al. 200^ coupled with laser guide star 
adaptive o ptics (AO; Wizinowich et al. 2006 van Dam 
et al.|20Q^ at the Keck 1 telescope on 2014 December 08 


Universal Time (UT). OSIRIS samples over 1000 spectra 
in a single observation, resulting in a three-dimensional 
data cube with a field of view of 0'.'32 x 1'.'28 spanning 
about 1600 spectral channels for the O'.'020/pixel plate 
scale with broad band filters. Our observations bene¬ 
fited from the recently upgraded “G3” grating ( jMieda 


et al. 2014). 

21X10441+2301 Aa {R = 14.2 mag) was used as a ref¬ 
erence tip tilt star for both binaries. The R-band mag¬ 
nitude of the laser guide star on the wavefront sensor 
remained steady at 9.3 mag during all the observations. 
Sky conditions were clear with excellent seeing between 
0.4-0.6" in the optical as recorded by the Differential Im¬ 
age Motion Monitor at the Canada-France-Hawaii Tele¬ 
scope. We obtained a total of 40 min of on-source in¬ 
tegration for 2M0441+2301 Aab in Jbb band (1.18-1.42 
/im), 30 min in Hbb band (1.47-1.80 /im), and 30 min 
in Kbb band (1.97-2.38 jam) over an airmass of 1.0 to 
1.2. For 2M0441+2301 Bab, we obtained 67 min in Jbb^ 
40 min in Hbb^ and 40 min in Kbb spanning airmasses 
between 1.0-1.6 (Table [B. The instrument rotator was 
set to 90° for 2M0441+2301 Aab and 120° for Bab, 
and the observations were executed in an ABBA pattern 
with 0.4-0.5" dithers along the long axis of the recon- 
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Fig. 2.— Near-infrared spectra of 2M0441-h2301 Aa, Ab, Ba, and Bb. Left: Flux-calibrated 1.2-2.4 p,m spectra. Major atomic and 
molecular features are labeled. The weak alkali lines such as K I at 1.25 p,m and triangular if-band (1.5-1.8 p,m) sha pes in the substellar 
components are signatures of low surface gravities and extreme youth. The near-infrared spectral types based on the|Allers & Liu| (|2013| 
classification scheme are M7 ± 1, M7 ± 1, and LI ± 1 for the Ab, Ba, and Bb components. All spectra have been smootlied trom 
to 17^1000 for clarity, normalized, and offset by a constant. The inset shows the unresolved 2MASS AT^-band image of the system together 
with our AO images of each subsystem in Ks with NIRC2. Contours for the Aab pair represent 50%, 20%, 7%, 4%, and 1.5% of the peak 
flux. Contours for the Bab pair represent 50%, 33%, 20%, 10%, and 5% of the peak flux. Note that the Aab image has been convolved 
with a Gaussian for better rendering. Right: The Pa/3 spectral region at the native (i7~3800) resolving power. 2M0441+2301 Aa and Ab 
show weak but significant emission with equivalent widths of -0.96 ± 0.02 A and -0.84 ± 0.11 A, respectively. The gray shaded region 
shows the 1 a flux density uncertainties. 


structed data cube. The AOV standard stars HD 19600 
and HD 35036 were targeted immediately before and af¬ 
ter both science objects with airmasses well matched to 
the associated science frames for each bandpass. 

Basic image reduction including bad pixel removal, flat 
fielding, pair-wise subtraction, wavelength calibration (in 
vacuum wavelengths), and assembly of the 2-D spectra 
into 3-D cubes was carried out using version 3.2 of the 
OSIRIS Data Reduction Pipeline with the most recent 
rectification matrices maintained by Keck Observatory. 
The current separation between 2M0441+2301Ba and Bb 
is 95 mas, which is only 1.8 diffraction limits at the IO¬ 
meter Keck telescope in K band. The fainter compan¬ 
ions in each binary pair (Ab and Bb) are easily visible in 
all median-collapsed data cubes, but their PSFs overlap 
with those of their primaries, so robust spectral extrac¬ 
tion free of systematics is not trivial. Our solution is to 
simultaneously fit analytical PSF models to each binary 
component by making use o f the nonlinear le ast-squares 
curve fitting routine MPFIT ( Markwardt|2009|) written in 
Interactive Data Language (IDL). We tested four models 
for each binary component: a single bivariate Gaussian, 
a single Lorentzian, the sum of a bivariate Gaussian and 
Lorentzian, and the sum of three bivariate Gaussians. 
These increasingly complex models are designed to mimic 
the extended PSF halo and a more compact diffraction- 
limited core. 

The PSFs of each component in our OSIRIS data have 
the same shape, remain at fixed relative positions in each 


cube, and both broaden in the same fashion at longer 
wavelengths, but their amplitudes vary with wavelength 
because of differing spectral slopes and absorption fea¬ 
tures. We therefore first fit each joint model to the bi¬ 
nary in the median-combined cube to establish the posi¬ 
tion of each model component, which is then held fixed 
for the PSF fit at each wavelength. All other parameters 
defining the PSF are allowed to vary, but the shape of 
each component (Aa and Ab, for example) is designed to 
remain mutually identical. Ultimately, we adopted the 
model composed of three Gaussians per binary compo¬ 
nent (i.e., six Gaussians per binary) because this pro¬ 
duced the best fit as measured by the lowest residual 
rms. Fourteen free parameters define the Gaussian am¬ 
plitudes, standard deviations, position angles, and over¬ 
all offset of the binary model (Figure and was used to 
fit our observations at each wavelength channel of every 
data cube and in all three filters for both 2M0441+2301 
Aab and Bab. The same model (except for a single com¬ 
ponent instead of a binary) was used to extract the spec¬ 
trum of the standard stars. 

The extracted spectra are scaled to a common level 
and the mean and standard error at each wavelength 
are adopted as the flux density and its associated un¬ 
certainty. Telluric correction was carried out with the 
xtellcor.general routine in the IRTF/Sp eX reduction 
package Spextool, which is written in IDL (|Vacca et al.| 
2003 Gushing et al.|2004 ). Each band was independently 


flux calibrated by hrst converting unresolved 2MASS 
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Fig. 3.— Spectral energy distributions for 2M0441+2301 Aab- 
Bab. Dark gray lines show the best matching BT-Settl mod¬ 
els (as indicated in the legend) to <3 /rm photometry assuming 
Ay =0.0 mag, which avoids disk emission at longer wavelengths. 
Shaded regions show the 1 cr and 2 cr uncertainties . Photometry 
origin ates from the literature (|Todorov et al.||2010| |Bulger et al.| 
[Mul l and our own NIRC2 and UblKlb observations Unresolved 
excess emission from one or more subdisks is present beyond about 
7 fim for the Aab component and 3 /rm for Bab, shown here as as¬ 
sociated with the primaries of each pair but which may also (or in¬ 
stead) originate from the companions. The flux-calibrated OSIRIS 
spectra are overplotted in black for comparison. 

photometry of each subsystem into the OSIRIS filter sys¬ 
tem by deriving photometric transformations as a func¬ 
tion of spectral type from synthetic photomet ry of MLT 


objec ts from the SpeX Prism Spectral Library (Burgasser 


2014). Binary flux ratios are measured during the mod¬ 


eling and spectral extraction of the science data and are 
used to calculate apparent magnitudes of each compo¬ 
nent in the OSIRIS filter system. Finally, a scale factor is 
derived to absolutely flux calibrate each bandpass of the 
spectrum to match the decomposed OSIRIS photometry 
(Figure [^. The signal-to-noise ratio per pixel is lower 
near the edges of each bandpass and is highest in the 
centers with median values across the entire 1.2-2.4/im 
spectrum of 112 for Aa, 45 for Ab, 41 for Ba, and 24 for 
Bb. 

In addition to our OSIRIS observations, we also ob¬ 
tained NIRC 2 images of 2M0441+2301 AabBab at the 
Keck II telescope on 2014 November 10 UT using Nat¬ 
ural Guide Star (NGS) AO. Sky conditions were clear 
with 0 .5-0.6" seeing in the optical throughout the obser¬ 
vations. We targeted the Aab component in K, J (on 
the Mauna Kea Observatory, or MKO, filter system), H 
(MKO), and Ks filters and the Bab component in H 
(MKO), and Ks filters with the narrow camera mode, 
resulting in a field of view of 10 " 2 xl 0 " 2 . A plate scale 
of 9.952 ± 0.002 mas pix“^ and sky orient ation on the 
detect or of 0.252 ± 0.009° are adopted from|Yelda et al. 


(2010). 2M0441-I-2301 Aa provided on-axis AO correc¬ 


tion for the Aab pair, and we offset by 12"3 for off-axis 
imaging of the Bab components. The images were dark 
subtracted, flat fielded, and cleaned of cosmic rays and 
bad pixels. Relative photometry and astrometry was 
measured by fitting three bivariate Gau ssians to each bi - 
nary component in each image following |Liu et al. (2010). 
We adopt the mean and standard deviation of these mea¬ 
surements, which are reported in Table 


3. RESULTS 

3.1. Spectral Classification 

The flux-calibrated spectra of 2M0441+2301 AabBab 
are shown in Figure The substellar components Ab, 
Ba, and Bb show strong atomic and molecular absorp¬ 
tion features characteristic of young late-M to early-L 
type brown dwarfs with low surface gravities. These in¬ 
clude shallow K I absorption lines at 1.25 jam and an 
angular spectral sha pe from 1.5 to 1.8 jam compared 
to old field objects (Allers & Liu 2013), which is pri¬ 
marily caused by weakened collision-induced absorption 
of molecular hydrogen at lower photos pheric pressures 
(Marley et al.|| 2012 | Barman et al.||2011|). 2M0441+2301 
Bb shows particularly pronounced features from water, 
carbon monoxide, iron hydride, and vanadium oxide. 
Index-based classifi cation of the substellar components 
(Allers & Liu||2013|) gives near-infrared spectral types of 
M7 ± 1, M7 ± 1, and LI ± 1 for 2M0441+2301 Ab, 
Ba, and Bb. We measure 1.2437 jam and 1.2529 jam K I 
equivalent widths of 3.5 ± 0.3 A and 0.0 ± 0.3 A for Ab; 
2.3 ± 0.3 A and 0.8 ± 0.5 A for Ba; and 3.4 ± 0.6 A 
and 5.5 ± 1.0 A for Bb, respectively. Aa and Ba have 
very low gravity (“VL-G”) classifications, while Bb has 
an intermediate-gravity classification (“INT-G”). 

3.2. Emission Lines and Accretion 

2M0441-1-2301 Aa and Ab show weak Pa^d emission in¬ 
dicating ongoing accretion from subdisks around both 
binary components. No significant Pa/3 emission is ev¬ 
ident in 2M0441+2301 Ba or Bb, and there is no Bry 
emission in any components. The emission line from the 
primary is superimposed on a slightly broader absorption 
troug h from unresolved Ti I, Fe I, and Ca I absorption 
lines ( Rayner et al.[2QQ9 ). We therefore define a pseudo¬ 
continuum level at 1.25 ± 0.01 x 10“^^ W m“^ /im“^ 
to compute the equivalent width of the emission line (- 
0.96 ± 0.02 A), line flux ( 1.20 ± 0.04 xlO”^^ W m-^), 
and line luminosity (log(I/Pa/ 3 /I/ 0 ) = -5.11 ± 0.09 dex). 
The corresponding accretion luminosity is log (Lace/A©) 
= -2.9 ± 1.0 dex using the empirical correlation b etween 
accretion lu minosity and Pa/d line luminosity from |Natta| 


et al.| (2004). We find a mass accretion rate of log(M) = - 
9.4 ± 1.1 M© yr“^ using M = 2 LaccR*/(GM^), where G 
is the gravitational constant and is the stellar mass. 
We adopt a stellar radius, R*, of 1.1 ± 0.2 R© based on 
the age and luminosity of the system and evolutionary 
models (Allard et al. 2012 ) and a stellar mass of 0.20^J;o 5 
M© (Section 3.4). 

For 2M0441±2301 Ab, we use the ± 0.01 /am region 
adjacent to the Pa/d line to define the pseudocontinuum 
level and measure a Pa/d equivalent width of-0.84 ± 0.11 
A. This corresponds to a line flux of 6.7 ± 0.9 x 10 “^^ 
W m“^, a line luminosity of log(Lpa/ 3 /L©) = -6.4 ±0.1 
dex, an accretion luminosity of log(Lacc/A©) = -4.7 ± 
1.3 dex, and a mass accretion rate of log(M) = -10.8 ± 
1.3 M(7) yr~ ^ assuming a radius of 0.41 ± 0.02 
et al.||2003 ) and a mass of 35 ± 5 Mj^p (Section 


Bar affe 


3.3. Atmospheric Model Fitting 

Atmospheric models can be used to infer the physical 
properties of the 2M0441±2301AabBab system such as 
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Fig. 4.— Comparison of 2M0441-h2301 AabBab with theoretical isochrones. Left: The BT-Settl models (|Baraffe et al.||2015|) spanning 
1, 2, and 3 Myr isochrones from right to left are consistent within 1 a with 2M0441+2301 AabBab in the Mh versus M—K near-infrared 
color-mag nitude diagrarn. All photometry has been converted to the MKO filter system. Comparison low-mass field objects are from the 
literature (|Dupuy & Liu|2012|), and symbols with dark borders represent y oung brown dwarfs an d giant planets as opposed to old stars and 
brown dwarts, l op tCigtit: llie 1 Myr isochrones from the Tucson models (|Burrows et al.|1997|) are generally consistent with the empirical 
isochrone traced by this quadruple system at the 2 cr level, though the Ba and Bb components are cooler than expected from these models 
by 200-300 K. 1 Myr, 5 Myr, 10 M yr, 100 Myr, and 1 G yr isochrones are plotted. Bottom Right: Same as the top right panel, but for the 
Lyon (Cond) evolutionary models (|Baraffe et al.||2003|). These models extend to 0.1 Mq and are consistent with 2M0441+2301 Ab, Ba, 
and Bb at the 2 a level. 1 Myr, 5 Myr, 10 Myr, foo IVlyr, and 1 Gyr isochrones are plotted. 


the components’ effective temperatures (Teff) and sur¬ 
face gravities (log g) by fitting synthetic spectra to our 
OSIRIS observatio ns. We first fit the BT-Settl models 
(lAllard et al.|2Q12|) w ith “CIFIST2011” solar abundances 
((Jaftau et al.||2UIUp to individual J, and K spectral 


regions and to the entire flux-calibrated 1.2-2.4 /im spec¬ 
trum of all components in the system. The grid spans 
1100 K to 4000 K in Teff(ATeff= 100 K) and 2.5 dex 
to 5.5 dex in log g (Alog = 0.5 dex [cgs units]). For 
each synthetic spectrum, we compute the scale factor 
that minimizes the value between the model and the 
data. This also gives the objects radius R since the scale 
factor is simply R^/dP^ where d is the distance to the 
object. We adopt the characteristic distance to Taurus 
of 145 ± 15 pc throughout this work. The models qual¬ 
itatively reproduce our near-infrared spectra with var¬ 
ious combinations of temperature and gravity, but the 
best-fit values to different bandpasses are generally not 
self-consistent and are therefore unreliable, yielding un¬ 
physical values that differ by up to 1300 K and 2.0 dex for 
the same object. This is likely a result of missing opac¬ 
ity sources, incomplete atomic and molecular line lists, 
and imperfect cloud prescriptions in the models. For 
example, although fits to the individual bandpasses of 


2M0441+2301 Bb provide good reduced chi-squared val¬ 
ues, the model effective temperatures range from 1600 K 
to 2500 K and the inferred radii span 1.6 to 5.4 Rjup- 
An alternative approach is to use photometry to con¬ 
strain the components’ effective temperatures since this 
spans a broader range of wavelengths and is less sensi¬ 
tive to the detailed physics and input line lists compared 
to our near-infrared spectra. After flux calibrating the 
models to the i^-band monochromatic flux density, we 
identify the best matches to the available photometry of 
all four objects spanning 0.4-3 gm. As with compar¬ 
isons to our spectra, the best-fitting models often have 
unphysical surface gravities. We therefore fix the sur¬ 
face gravitie s to the predicted values from the |Baraffe| 


et al. (2003) evolutionary models because gravity is only 
minimally influenced by model assumptions. We adopt 
values of log = 4.0 dex for Aa and log ^ = 3.5 dex 
for Ab, Ba, and Bb and compute values for a range 
of effective temperatures (1500 K to 4000 K). The best 
matches have effective temperatures of 3400 K for Aa, 
2800 K for Ab, 2100 K for Ba, and 1800 K for Bb. 

The results for Ba and Bb are up to 500 K cooler 
than recent spectral type -effective temperature calibra - 
tions for young stars from Herczeg & Hillenbrand (2014), 
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which are derived by fitting the same grid of BT-Settl 
atmospheric models we use in this work to optical spec¬ 
tra of T Tauri stars. This discrepancy may reflect the 
paucity of photometric observations for 2M0441+2301 
Ab, Ba, and Bb shortward of J band where the SED 
turns over and provides the best constraint on Teff. The 
reduced values (y^) are much larger than unity 
^ 10 to 300), implying the photometric uncertainties are 
underestimated and/or there are systematic errors in the 
models. Since these data do not appear to be standard 
deviates of the models, maximum likelihood methods for 
establishing parameter confidence regions like constant 
boundaries about the minimum values are not 
valid here. We therefore adopt a single grid step (100 K) 
as an estimate of the uncertainty in Teff for Aa and two 
grid steps (200 K) for the substellar objects. Note that 
we have assumed zero reddening for these fits, but red¬ 
dening up to Ay = 0.5 mag does not change the fitting 
results by more than 1 a. Additional resolved optical 
photometry would provide better constraints on Teff and 
log which in turn would enable more rigorous tests of 
substellar evolutionary models. 


3.4. Luminosities and Masses 

Integrating the models and adopting the mean Taurus 
distance gives luminosities of \og{L/LQ) = -0.86 ± 0.09 
dex, -2.13 ± 0.09 dex, -2.46 ± 0.09 dex, and -3.03 ± 
0.09 dex for Aa, Ab, Ba, and Bb, respectively. Assuming 
a uniform likelihood of ages between 1 to 3 Myr, the cor¬ 
responding masses of the substellar components are 35 ± 
5 Mjup for Ab, 19 ± 3 Mj^p for Ba, and 9.8 ±1.8 Mj^p 
for Bb from interpo lated substellar evolutionary models 
(Baraffe et al.||2003). We adopt a mass of 0.20lo 05 Mq 
for the primary star Aa based on recent low mass stel¬ 
lar evolutionary models ( Baraffe et al.||2Ql^. The mass 
we fin d for Bb is somewhat greater than |Todorov et al.| 
( 2010 ) owing to our slightly higher luminosity measure¬ 
ment and broader assumed age range of 1 to 3 Myr. 

Altogether this confirms that 2M0441±2301 Bb is sim¬ 
ilar in mass and separation to directly imaged planets like 
HR 8799 bc de dMarois et ah' 2008| Marois et al. 2010) 


and P Pic b (Lagrange et al. 


2010| ). 'Lhe mass ratios of 


these systems, however, ditter by nearly two orders of 
magnitude; the HR 8799 and /3 Pic systems have mass 
ratios of about 150-220, but for the 2M0441±2301 Bab 
binary it is only about 1.9. Together with environmen¬ 
tal clues such as apparently coplanar multiple planets 
and an unusually massive and extended debris disk, the 
high mass ratios of HR 8799 and P Pic point to an al¬ 
ternative formation route compared with 2M0441±2301 
Bab (see Section 0. Finally, these results also confirm 
that 2M0441±230T AabBab is the lowest-mass quadru¬ 
ple system currently known, with a total mass of about 
0.26 Mq. 

4. DISCUSSION 

As a coeval and presumably chemical homogeneous 
quadruple system spanning the stellar to planetary 
mass regimes, 2M0441±2301 AabBab also provides a 
unique opportunity to test pre-main sequence and gi¬ 
ant planet evolutionary models similar to other high- 
order stell ar/substellar multipl e systems in Taurus like 


GG Tan (White et al 
location of 21VIO441±230 


19991 

Tfa 


Figure shows the 
oBab in color-magnitude 


and Hertzsprung-Russell diagrams. Evolutionary model 
isochrones spanning 1 to 3 Myr are broadly consistent 
with the data at the 1 to 2 a levels. 

We detect orbital motion for the first time in both the 
Aab and Bab binaries. 2M0441±2301Aab has undergone 
significant changes in separation of 1.4 ± 0.2 mas yr“^ 

m onrrlo A ^ ^ O 09° Trr--! 


2008 (Kraus et al. 2011 

Todorov et al. 

2010 

Todorov 

et al.|2U14l). 2MU441+230 

11 Bab is moving by -] 

L.5 ± U.5 


Innately, the long orbital periods of about 300-400 years 
for each subsystem means that dynamical mass measure¬ 
ments will not be possible for at least several decades. 

The wide separation between 2M0441±2301 Aab and 
Bab, the hierarchical orbital architecture of the system, 
and the low mass ratios of the Ab, Ba, and Bb com¬ 
ponents relative to the primary star Aa argue against 
formation by core accretion or stellar capture. Simula¬ 
tions of fragmenting disks can result in close substellar 
companions and ejected brown dwarf-planetary mass ob¬ 
ject binaries, but such low mass quadruple systems with 
a near ejection of one pair is difficult to reproduce with 
this formation channel (Li et al.||2015 ). Instead, the 2±2 
hierarchy of 2M0441±23(J1 AabBab resembles the con¬ 
figuration of the most common stellar quadruples, which 
outnumber 3±1 hierarchies in planetary-like configura- 
tio ns by at least 4 t o 1 among Sun-like stars within 25 
pc (Tokovinin 2014). High order multiple systems like 
2M0441±23U1 AabBab provide important clues about 
the masses, critical densities, characteristic sizes, and 
opacity limits of fragmenting molecular cloud cores. This 
is particularly true at young ages and in low density star 
forming regions before significant dynamical unfolding 
has occurred and where interactions with nearby stars 


are m inimal (Reipurth & Mikkola 2015 Pineda et al 


2015). Although the details and universality of quadru¬ 


ple star formation remain unclear, one possible formation 
route for this system is through two primary fragmenta¬ 
tion events which each subsequently underwent a sec¬ 
ondary fragmentation following H 2 dissociation, result¬ 
ing in two short-pe riod binaries on weakly bound or bits 
about each other (Whitworth & Stamatellos 2006). If 
the current separation of at least 1800 AU reflects the 
nascent architecture of the system, then the binary pairs 
have undergone fewer than 20 full orbits around each 
other in 3 Myr. 

The lower limit of cloud fragmentation has been ex¬ 
tensively explored with increasingly sophisticated sim¬ 


ulatio ns of collapsing turbulent molecular clouds (Bate 


2009). Opacity-limited fragments between about 3 to 


10 Mjup can be produced if accretion is rapidly halted 
through dyn amical encounters or ejecti ons out of dense 
cloud cores (Padoan & Nordlund 2004). This is gener¬ 
ally consistent with empirical determinations of the ini¬ 
tial mass function, which show a smooth and continuous 
distribution down to 5 to 10 Mjup and sugge st a com¬ 


mon origin with most stars and brown dwarfs (Ghabrier 


2003). This mass range also coincides with the lowest- 


mass companions ima ged around stars and brown dwarfs 
(Brandt et al. 2014). However, the formation path- 
way of individual systems is often ambiguous since free- 
floating planetary-mass objects could represent ejections 
from planet-planet scattering events, and directly imaged 
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planets could be scattered or in situ products of disk in¬ 
stability or core accretion. Confirming the low mass of 
2M0441+2301 Bb demonstrates that fragmenting molec¬ 
ular clouds can create companions which overlap in mass 
with bona fide planets formed in disks, implying that 
some of the directly imaged planets orbiting stars and 
brown dwarfs may represent the tail end of turbulent 
fragmentation rather than the tip of the iceberg of planet 
formation. 
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and SIMBAD database operated at CDS, Strasbourg, 
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tral Libraries, maintained by Adam Burgasser at 
http: / / pono.ucsd.edu / ^adam/browndwarfs / spexprism, 
Finally, mahalo nui loa to the kama‘aina of Hawaifi for 
their support of Keck and the Manna Kea observatories. 
We are grateful to conduct observations from this 
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TABLE 1 

Relative Photometry and Astrometry 


Name 

UT Date 

Inst. 

Ax Exp. 

Filt. 

Sep. 

RA. 

A mag 

rriA 

rriB 

FWHM 




Time (s) 


(") 

(°) 




(mag) 

mag 

(mas) 

2M0441Aab 

2014 Nov 10 

NIRC2 

7 X 15 

Y 

230 ± 2 

79.3 ± 0.6 

3.16 


0.16 

11.37 ± 0.02 

14.53 

± 0.15 

96 ± 29 

2M0441Aab 

2014 Nov 10 

NIRC2 

11 X 5 

J 

234 ± 2 

79.4 ± 0.6 

2.96 


0.06 

10.77 ± 0.04 

13.73 

± 0.07 

65 ± 12 

2M0441Aab 

2014 Nov 10 

NIRC2 

7x5 

H 

233.0 ± 1.6 

79.9 ± 0.2 

3.11 


0.07 

10.19 ± 0.03 

13.30 

± 0.07 

65 ± 17 

2M0441Aab 

2014 Nov 10 

NIRC2 

6x5 

Ks 

232.0 ± 1.3 

79.1 ± 0.3 

2.77 


0.13 

9.94 ± 0.02 

12.71 

± 0.12 

64 ± 5 

2M0441Aab 

2014 Dec 08 

OSIRIS 

8 X 300 

Jbb 



3.3 


0.3 

10.70 ± 0.05 

14.0 

± 0.3 

54 ± 5 

2M0441Aab 

2014 Dec 08 

OSIRIS 

6 X 300 

Hbb 



3.14 


0.18 

10.20 ± 0.03 

13.34 

± 0.17 

58 ± 3 

2M044lAab 

2014 Dec 08 

OSIRIS 

6 X 300 

Kbb 



2.81 


0.10 

9.96 ± 0.02 

12.77 

± 0.10 

54 ± 1 

2M0441Bab 

2014 Nov 10 

NIRC2 

7 X 60 

H 

93 ± 2 

132.6 ± 1.1 

1.45 


0.14 

14.03 ± 0.05 

15.48 

± 0.12 

88 ± 10 

2M0441Bab 

2014 Nov 10 

NIRC2 

7 X 30 

Ks 

97.5 ± 1.0 

132.0 ± 0.6 

1.24 


0.05 

13.46 ± 0.03 

14.70 

± 0.05 

65 ± 4 

2M0441Bab 

2014 Dec 08 

OSIRIS 

10 X 400 

Jbb 



1.7 


0.2 

14.57 ± 0.06 

16.26 

± 0.17 

67 ± 7 

2M0441Bab 

2014 Dec 08 

OSIRIS 

8 X 300 

Hbb 



1.42 


0.19 

14.06 ± 0.05 

15.47 

± 0.15 

88 ± 9 

2M0441Bab 

2014 Dec 08 

OSIRIS 

8 X 300 

Kbb 



1.35 


0.17 

13.50 ± 0.05 

14.85 

± 0.14 

79 ± 10 



